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for  the  rapid  evaluation  of  the  potential  and  force  fields  in  systems  of  particles  whose  interactions 
are  Coulombic  or  gravitational  in  nature.  For  a  system  of  N  particles,  an  amount  of  work  of 
the  order  0(N *)  has  traditionally  been  required  to  evaluate  all  pairwise  interactions,  unless  some 
approximation  or  truncation  method  is  used.  The  algorithm  presented  here  requires  an  amount  of 
work  proportional  to  N  to  evaluate  all  interactions  to  within  roundoff  error,  making  it  considerably 
more  practical  for  large-scale  problems  encountered  in  plasma  physics,  fluid  dynamics,  molecular 
dynamics  and  celestial  mechanics.  _ . 
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1.  Introduction 


The  evaluation  of  Coulombic  and  gravitational  interactions  in  large-scale  ensembles  of  particles 
is  an  integral  part  of  the  numerical  simulation  of  a  large  number  of  physical  processes.  Typical 
examples  include  celestial  mechanics,  plasma  simulations,  the  vortex  method  in  fluid  dynamics, 
and  the  solution  of  the  Laplace  equation  via  potential  theory  (see  [10],  [5],  [14],  [2]).  In  such  cases, 
the  potential  has  the  form 


®  ^ external  "b  local  "b  ^/ar>  (1-1) 

where  $  local  is  a  rapidly  decaying  function  of  distance  (such  as  the  Van  der  Waals  potential  in 
chemical  physics),  $  external  is  a  function  which  is  independent  of  the  number  and  relative  positions 
of  the  particles  (such  as  an  external  gravitational  field),  and  $/ar  is  Coulombic  or  gravitational. 

In  the  numerical  evaluation  of  fields  of  the  form  (1.1),  the  cost  of  computing  the  terms  $ external 
and  $i0cai  is  of  the  order  O(N),  where  N  is  the  number  of  particles  in  the  ensemble.  Indeed,  $ external 
is  evaluated  separately  for  each  particle,  and  $!<*<,/  decays  rapidly,  involving  the  interactions  of  each 
particle  with  a  small  number  of  nearest  neighbors.  Unfortunately,  evaluation  of  the  term  $/or,  if 
done  directly,  requires  order  0(iV2)  operations,  since  the  Coulombic  potential  decays  slowly,  and 
the  interactions  between  each  pair  of  particles  have  to  be  taken  into  account.  In  many  situations, 
in  order  to  be  of  physical  interest,  the  simulation  has  to  involve  thousands  of  particles  (or  more), 
making  the  estimate  0(N2)  excessive  in  some  cases,  and  prohibitive  in  others. 

Many  three-dimensional  processes  are  more  or  less  adequately  described  by  two-dimensional 
models,  and  this  fact  is  widely  used  in  computer  simulations.  From  the  computational  point  of  view, 
reduction  of  the  dimensionality  of  the  problem  has  several  major  advantages:  fewer  particles  are 
normally  required  to  obtain  a  physically  meaningful  model  of  a  two-dimensional  process  than  that 
of  its  three-dimensional  counterpart,  the  numerical  techniques  for  calculations  in  two  dimensions 
are  better  developed  and  easier  to  implement,  and  finally,  the  display  and  interpretation  of  three- 
dimensional  results  pose  problems  almost  non-existent  in  two  dimensions.  On  the  other  hand, 
certain  processes  in  the  physical  world  simply  can  not  be  approximated  by  two-dimensional  models. 
In  such  cases,  full  three-dimensional  simulations  have  to  be  performed,  with  the  help  of  appropriate 
numerical  tools. 

We  restrict  our  attention  now  to  the  evaluation  of  Coulombic  or  gravitational  potentials,  and 
begin  by  observing  that  multipole  expansions  have  been  a  popular  tool  in  the  analysis  of  such 
problems  for  more  than  100  years  ([11], [12]).  They  are  routinely  used  to  represent  the  fields  due  to 
collections  of  charges  (masses)  in  regions  of  space  removed  from  the  source  positions.  However,  these 
classical  techniques  are  generally  inapplicable  in  situations  where  the  charges  and  the  points  where 
the  fields  are  to  be  evaluated  are  not  separated.  The  use  of  multipole  expansions  for  the  evaluation 
of  fields  in  such  a  situation  was  reported  in  [14],  where  the  solution  of  the  Laplace  equation  by 
means  of  boundary  integrals  is  discussed.  The  method  described  requires  the  evaluation  of  the 
potential  field  induced  by  a  collection  of  N  charges  lying  on  a  curve  at  each  of  the  charge  locations, 
and  a  multipole- based  algorithm  is  used  to  carry  out  this  calculation  in  order  O(N)  operations.  In 
[7], [8],  we  introduced  a  generalization  of  this  method  for  the  rapid  evaluation  of  the  potentials  and 
forces  in  large-scale  two-dimensional  systems  of  particles  randomly  distributed  in  a  square  domain. 
This  fast  multipole  algorithm  requires  an  amount  of  work  proportional  to  N  to  evaluate  to  within 
round-off  error  all  pairwise  interactions  in  a  system  of  iV  charges.  In  [4],  an  adaptive  version  of 
the  algorithm  of  [7], [8]  was  introduced,  whose  CPU  time  requirements  are  proportional  to  JV  and 
independent  of  the  statistics  of  the  charge  distribution. 

This  paper  reports  the  theoretical  apparatus  underlying  the  three-dimensional  version  of  the 
fast  multipole  algorithm.  For  a  full  description  of  the  method  in  three  dimensions,  as  well  as  a 


detailed  discussion  of  numerical  results  produced  by  implementations  of  the  method  in  two  dimen¬ 
sions,  we  refer  the  reader  to  [13].  Other  aspects  of  the  problem  discussed  in  [13]  include  applications 
of  the  method  to  several  problems  in  physics,  chemistry,  and  biology. 

While  in  [7], [8]  the  analytical  foundation  is  the  theory  of  holomorphic  functions,  here  it  is  the 
theory  of  spherical  harmonics.  Conceptually,  the  transition  is  quite  straightforward.  However,  a 
number  of  technical  problems  had  to  be  overcome  before  an  efficient  three-dimensional  algorithm 
could  be  constructed.  In  particular,  we  formulate  and  prove  two  generalizations  of  the  classical 
addition  theorem  for  the  Legendre  polynomials  (Theorems  2.3  and  2.4  of  this  paper)  that  appear 
to  have  been  previously  unknown. 

Remark:  Several  other  approaches  have  been  used  to  reduce  the  cost  of  the  evaluation  of  Coulom- 
bic  potential  fields.  For  a  detailed  discussion  of  these  algorithms,  we  refer  the  reader  to  [7],  and  to 
the  original  works  [10],  [l],  [2]. 


2.  Physical  and  Mathematical  Preliminaries 


In  this  paper,  we  will  be  considering  three-dimensional  physical  models,  consisting  of  a  set 
of  jV  charged  particles  with  the  potential  and  force  obtained  as  the  sum  of  pairwise  interactions 
via  Coulomb’s  law.  If  a  point  charge  of  unit  strength  is  located  at  the  origin,  then  for  any  point 
P  =  (x,  y,  z )  €  R3  with  ||.P||  =  r  ^  0,  the  potential  and  electrostatic  field  due  to  this  charge  are 
described  by  the  expressions 

$  =  -  and 
r 


respectively. 

Using  spherical  coordinates,  suppose  now  that  the  unit  charge  is  located  at  a  point  Q  — 
( p,a,0 ),  not  the  origin.  The  potential  at  a  point  P  Q  is,  of  course,  the  inverse  of  the  distance 
117*  -  Q\\  —  r1  (Figure  1).  There  is  a  well-known  series  expansion  for  the  potential  at  P  =  ( r,9,4> ), 
in  terms  of  its  distance  from  the  origin  r.  Letting  qr  be  the  angle  between  the  vectors  P  and  Q ,  we 
have 


(2.1) 


where  u  =  cosqi  and  P„(u)  is  the  Legendre  polynomial  of  degree  n.  Equation  (2.1)  is  generally 
referred  to  as  a  multipole  expansion,  and  describes  the  far  field  due  to  a  charge  at  Q ,  since  the 
condition  for  its  validity  is  that  r  >  p. 

There  is  a  duality  inherent  in  the  situation  depicted  in  Figure  1,  namely  that  if  the  locations 
of  the  charge  (Q)  and  the  evaluation  point  (P)  were  interchanged,  then  the  field  at  P  would  still 
be  described  by  p.  In  this  case,  so  long  as  r  <  p,  we  may  write 


Equation  (2.2)  is  valid  only  in  the  open  sphere  centered  at  the  origin  with  radius  p.  and  we  will 
refer  to  such  a  description  of  the  potential  field  as  a  local  expansion. 

The  following  observation  pertaining  to  Legendre  polynomials  will  be  needed  below.  Its  proof 
can  be  found  in  most  standard  textbooks  (see,  for  example,  [11]). 
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Lemma  2.1.  Let  u  €  R,  with  |u|  <  1.  Then  Pn(u)  <  1. 

From  this  fact,  it  is  straightforward  to  obtain  the  following  two  error  bounds. 

Lemma  2.2.  Suppose  that  a  charge  of  strength  q  is  located  at  the  point  Q  =  (p,  a,/?),  and  that 
P  =  (r,9,<p)  €  R3,  with  ||P  —  Q\\  =  r1  and  r  >  p.  Letting  7  be  the  angle  between  the  two  points , 
we  have  an  error  bound  for  the  multipole  expansion  (2.1)  of  the  form 


<-vh(T'  ■ 


n=0 


(2.3) 


Similarly,  when  r  <  p,  we  have  an  error  bound  for  the  local  expansion  (2.2)  of  the  form 


1  -V'  q  r" 

r 1  2-*  pn+1 
n=0 


Pn(cos'y) 


p-  r  \p, 


,  p+i 


(2.4) 


Examining  equation  (2.1)  yet  again,  we  note  that  the  resulting  series  depends  on  the  relative 
coordinates  of  the  two  particles.  If  another  such  series  were  developed  for  a  second  source  at  the 
point  Q1,  they  would  have  to  be  evaluated  independently.  We  will  therefore  need  a  more  general 
approach  to  the  solution  of  potential  problems,  which  will  allow  us  to  develop  asymptotic  expansions 
of  the  field  due  to  arbitrary  sets  of  charges. 


2.1.  Spherical  Harmonics 

The  development  of  a  general  expansion  describing  potential  fields  in  three  dimensions  is 
carried  out  by  considering  the  Laplace  equation  itself,  which  characterizes  the  behavior  of  such 
fields  in  regions  of  free  space.  Using  spherical  coordinates,  the  Laplace  equation  takes  the  form 

l_d_  (  2dd>\  1  d  (  .  ££\  1  d2$  _ 

r2  dr  \  dr  )  +  r2  sin  9  dO  \  m  dO  )  r2  sin2  0  d<t> 2 

The  standard  solution  of  this  equation  by  separation  of  variables  results  in  an  expression  for  the 
field  as  a  series,  the  terms  of  which  are  known  as  spherical  harmonics. 

00  n  /  \ 

*(r,M)  =  £  £  (c-^  +  ^J-W,*)  (2.5) 

n=0  m=-n  '  ' 

In  the  above  expansion,  the  terms  Y™{0,  4>)rn  are  usually  referred  to  as  spherical  harmonics  of  degree 
n,  the  terms  Y™{9,  <t>)/rn+l  are  called  spherical  harmonics  of  degree  —  n  -  1,  and  the  coefficients 
LJJ*  and  Af™  are  known  as  the  moments  of  the  expansion. 

Remark:  It  is  obvious  that  in  a  far  field  (multipole)  expansion,  the  coefficients  L\ ”  must  be  set  to 
zero  in  order  to  satisfy  the  condition  that  the  field  decay  at  infinity.  In  a  local  expansion  (which  is 
to  be  analytic  in  a  sphere  centered  at  the  origin),  it  is  clearly  the  coefficients  XI™  which  must  be 
set  to  zero. 

Lemmas  2.3  -  2.5  below  are  well-known.  Their  proofs  can  be  found,  for  example,  in  [9]  or  [16). 

Lemma  2.3. 

*;°(M)  _  4o  ^_(i\ 

rn+l  dzn\r)  ' 
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For  m  >  0,  we  have 


mc4 


W'Q)  =  A"'.iJL  +  ,ir  ( AV“m  ( i\ 

rn+ 1  n  dx  dy  \dz)  \r ) 

CTj)  _  4m .  ( A  _  ,±v»  (i\ 

rn+l  ‘  "  3x  dy  \dz  )  \r ) 


where 


-*C  = 


(-1)" 


n  -  m  !  n  +  m  ! 


Remark:  The  standard  definition  of  the  functions  yrJn(^,<^)  includes  a  normalization  factor  of 

\j2n  4-  1/4tt.  (2.7) 

We  will  consistently  use  the  definition  given  above.  That  is,  the  coefficient  (2.7)  will  always  be 
omitted. 

Since  certain  differential  operators  arise  frequently  in  discussions  of  spherical  harmonics,  we 
introduce  the  following  notation. 

Definition  2.1.  The  operators  .  d and  d,  are  defined  by  the  expressions 

d  d 

d±  —  —  ±  i  —  and 
dx  dy 

Lemma  2  4  It  ;  i>  a  harmonic  function,  then 

d+d-  [<t>)  =  -d]{<{>) 


Lemma  2  5  f  r  im  n  s  rn  >  0. 

'  W"  (  \  )  =  (-!)"(«  -  H)!^  PjT'(coS9)  ■  e±im*  , 

where  the  term  /^("(cos#)  denotes  the  associated  Legendre  function  of  degree  n  and  order  m. 
Combining  lemmas  2.5  and  2.3.  we  have  a  useful  expression  for  the  spherical  harmonics: 

V„m(^)  =  Y  ~  j~jjj  ^m|(cos0)e'^.  (2.8) 

2.2.  The  Field  Due  to  Arbitrary  Distributions  of  Charge 

We  will  need  a  well-known  result  from  the  theory  of  spherical  harmonics,  which  is  usually 
referred  to  as  the  Addition  Theorem. 


•41 


I 

I 

i 

I  Theorem  2.1.  (Addition  Theorem  for  Legendre  Polynomials )  Let  P  and  Q  be  points  with  spherical 

coordinates  (r,6,<p)  and  (p,a,/?),  respectively,  and  let  ~t  be  the  angle  subtended  between  them 
(Figure  1).  Then 

n 

Pn( cos7)=  E 

m=— n 

It  is  a  straightforward  matter  now  to  form  a  multipole  expansion  describing  the  far  field  due 
to  a  collection  of  particles. 

Theorem  2.2.  (Multipole  Expansion).  Suppose  that  m  charges  of  strengths  {g,-,  :  =  l,...,m}  are 
located  at  the  points  {Q,  =  (p,-,  a,-,  /?;),  i  =  1, ...,  m},  with  |p,|  <  a.  Then  for  any  P  =  (r,  9,  <t>)  €  R3 
with  r  >  a,  the  potential  4>{P)  is  given  by 

HP)  =  E  E  (2.9) 


where 

m 

A/„m  =  E?'  P? 

1=1 

Furthermore,  for  any  p  >  1, 


W  -  E  E 


AC 


< 


(;)'+1  ■ 


(2.10) 


(2.11) 


where 

A  =  fki|.  (2.12) 

i=i 

Proof.  Let  us  first  consider  the  contribution  from  a  single  charge  g,  located  at  (p,-,  a,-,  /?,■).  From 
formula  (2.1)  and  the  Addition  Theorem  for  Legendre  Polynomials,  we  have 

H  =  E  7 '  ■P"(cosT') 

n=0 

o°  n  r  n 

=  E  E 

n=0 m=-n 

The  coefficients  JVC  in  equation  (2.10)  are  then  obtained  by  superposition.  The  error  bound  is 
an  immediate  consequence  of  (2.3),  the  triangle  inequality,  and  the  fact  that  the  ratios  p,/r  are 
bounded  from  above  by  a/r. 

I 

Before  proceeding  with  the  further  development  of  the  theory  of  spherical  harmonics,  we 
will  demonstrate  with  a  simple  example  how  multipole  expansions  can  be  used  to  reduce  the 
computational  complexity  of  the  evaluation  of  potential  fields.  Suppose  that  a  collection  of  m  point 
charges  of  strengths  {g,-,  i  =  l,...,m}  are  located  at  the  points  {Q,  =  (p,-, a,-, /?,-),  i  =  l.-  -.m), 
and  that  {Pj  =  ( ry,#;- ,  4>j),  j  =  1,  •  •  ,  n}  is  another  set  of  points  in  R3  (Figure  2).  We  say  that  the 


rn+ 1 


Y?(9,<t>) 


5 


sets  { Qi }  and  {Pj}  are  well-separated  if  there  exist  points  Po,Qo  €  R3  and  a  real  number  a  >  0 
such  that 

HQ;  -  Qo||  <  a  for  i  =  , 

|| Pj  -  Poll  <  a  for  j  =  1,  n  ,  and 
||Qo  -  Po||  >  3a. 

In  order  to  obtain  the  potential  at  each  of  the  points  Pj  due  to  the  charges  at  the  points  Q,  directly, 
we  could  compute 

m 

X^<?.(Pi)  for  j  =  l . n.  (2.13) 


This  requires  order  n-  m  work  (evaluating  m  fields  at  n  points).  Suppose,  on  the  other  hand,  that  we 
first  compute  the  coefficients  of  a  p'h-degree  multipole  expansion  of  the  potential  due  to  the  charges 
...,9m  about  Qo,  using  Theorem  2.2.  This  requires  a  number  of  operations  proportional  to 
m  ■  p 2.  Evaluating  the  resulting  multipole  expansion  at  all  points  Pj  requires  order  n  p2  work,  and 
the  total  amount  of  computation  is  of  the  order  0(m  ■  p2  +  n  ■  p2).  Moreover,  by  (2.11), 


Z>Q.(P>)  £  5Z  ||p  _Q0||n+i 

i=l  n=Om=-n  11  1  ^  " 


Y?{0h4>i)  <  Q)’ 


and  in  order  to  obtain  a  relative  precision  e  (with  respect  to  the  total  charge),  p  must  be  of  the 
order  —  log2(e).  Once  the  precision  is  specified,  the  amount  of  computation  has  been  reduced  to 

O(m)  +  O(n)  , 

which  is  a  significant  reduction  in  complexity  when  compared  with  the  direct  method. 


2.3.  Translation  Operators  and  Error  Bounds 

As  in  the  two-dimensional  case,  the  principal  analytical  tools  required  by  the  fast  algorithm  are 
certain  translation  operators,  acting  on  both  multipole  and  local  expansions.  In  order  to  develop 
the  necessary  formulae  for  these  procedures,  we  will  need  the  following  three  theorems,  which  can 
be  viewed  as  generalizations  of  the  classical  addition  theorem  for  Legendre  polynomials.  While  a 
somewhat  different  form  of  Theorem  2.5  below  can  be  found  in  the  literature  ([6],[3],[l5]),  Theorems 
2.3  and  2.4  appear  to  be  new.  The  following  theorem  describes  a  formula  for  the  expansion  of  a 
spherical  harmonic  of  negative  degree  about  a  shifted  origin. 

Theorem  2.3.  (First  Addition  Theorem)  Let  Q  =  (p,a,0)  be  the  center  of  expansion  of  an 
arbitrary  spherical  harmonic  of  negative  degree.  Let  the  point  P  =  (r, &,</>),  with  r  >  p.  and 
P-  Q  =  (P,0', </>').  Then 


n'+l 


=  Y'  J%'  A?  A$  Pn  Yn~m (a, (3)  Y™*?' 

/  -■  /  >  *m+m'  rn+n'+l 

An+n' 


n= 0 m=-n 


where 


Mi: 


_l)m*"<lm'l,M),  ifm  m!  <  0; 

,  otherwise. 


(2.14) 


Proof.  Making  use  of  equation  (2.1),  the  Addition  Theorem  for  Legendre  polynomials,  and  Lemma 
2.3,  we  observe  that 
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UP-QII  r1 


\ - 


=EE^rMi^ 

n=0  m=  — n 


=  E(  £  Pn  ■  ^"W)  ■  4T  '  alm|arw  (;)  + 

n=0  m=— /i  v  ^ 

£  P"  •  Y~m{a,0)  ■  ■  d%d?-m  (i)  )  . 


We  now  consider  three  separate  cases. 
Case  I  :  m'  —  0.  From  Lemma  2.3, 


YZ(8',4>')  _  .0  d„'(l 

rin'+l  An'  l  jJ 


(*)• 


Combining  (2.15)  and  (2.16),  we  obtain 


Y”fn'lt']  «  £(  £  ^  ^Lm|5rn'-|m|  (;)  ■ 

n=0  m=— n  '  7 

E  o"  I'.-fa./J)  A°„,  •  at  W“'-"  (i)  ) 


_V-  A  t P"  Y„-"{a,l3)  A°,  A?\  Y^(e,t) 

^  ^  \  A™,  /  r"+"'+1  ’ 

n=0  m=— n  \  n+n  / 


(2.15) 


(2.16) 


where  the  last  equality  is  obtained  by  another  application  of  Lemma  2.3. 

Case  II  :  m'  <  0.  Using  Lemma  2.3  again, 

r/n'+l  ”  ~An'  d-  *.  { r ,) 

=  £(  £  Pn  ■  Ynm{<*,0)  ■  4S'  ■  K  ■  dLm'l+Ha"+"'-H-|m'l  + 

n=0  m=-n 

E  (>”  •  Yn-"(a,0)  ■  .45’  ■  .45  ■  aKiaja;-*-”'-"-1’"'1  ( ' )  ) 


=  ££ 
n=0  m=— n 


^  / JS'  ■  45’  ■  45  ■  ■  r„— (°,gA  jgT  IM) 


4  m+m' 
An+n' 


_n+n'+l 


where 


/  1,  if  m  <  0; 

.  -  |  (-1)>»<'i(|mW,  if  m  >  0. 


To  obtain  the  last  equality,  for  the  terms  with  m  >  0,  we  have  used  Lemma  2.4  to  annihilate 
whichever  of  the  operators  d_  and  d+  occurs  less  frequently. 


Case  III  :  m'  >  0.  From  Lemma  2.3, 


~  n'+l  “An'  c'+c'* 

=  £(  H  /»"  '  Vn"m(a'^)  4?'  •  a“'aLm|a”+n'~|m|_m'  (;)  + 

n=0  m=— n  '  ' 


where 


£  p"  •  •  A-  ■  d™+m'd1 


1  an+n'-m-m'  /  ^ 


=  E  E 

n=0  m=-n 


J£' ■  A#  •  A”  ■  p"  Yn-m(a,l3)\  Y^f{e^) 


Am+m' 
n+n' 


rn+n'+ 1  ’ 


,/  _  f  1,  if  m  >  0; 

1  _  \  (_i)m«"("‘,.M),  if  m  <  0. 


As  before,  for  the  terms  with  m  <  0,  we  have  used  Lemma  2.4  to  annihilate  whichever  of  the 
operators  and  d+  occurs  less  frequently. 


The  second  addition  theorem  yields  a  formula  for  converting  a  spherical  harmonic  of  negative 
degree  (a  multipole  term)  with  respect  to  one  origin  into  a  local  expansion  about  a  shifted  origin. 

Theorem  2.4.  (Second  Addition  Theoremj  Let  Q  —  (p,a,0)  be  the  center  of  expansion  of  an 
arbitrary  spherical  harmonic  of  negative  degree.  Let  the  point  P  =  (r,9,<fr),  with  r  <  p  and 
P-  Q  =  {r’,d',<t>').  Then 


Y$\d\<t>')  _  f,  ^  -C'  -  Ay  ■  A™!  ■ 

pn'  +  l  Z-t  2-s  „n+n'+l  .  /tm-m' 

n=0m=-n  “  ^n+n' 


Ynm(e.4>)r "  , 


where 


r's 

*  u-ir, 


(—  l)m*n(|m,|,|m|)^  ifm  m>  >  Q. 

,  otherwise. 


(2.17) 


Proof.  We  first  let  (xp,yp,  zp)  and  (x<j,  Vq,  zq)  denote  the  Cartesian  coordinates  of  the  points  P 
and  Q,  respectively.  Then 


dxp  yr1 )  dxQ  \r' J 
dyP  \r' J  dyQ  \r' J 

A-(L)=-±_(L)  . 

dzp  \r' )  dzQ  \r' J 


(2.18) 


We  will  denote  by  d+p ,  d-p,  dZp,  d+Q,  d-Q ,  dZQ,  the  differential  operators  given  by  Definition  2.1. 
with  respect  to  the  indicated  variable  point. 
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Combining  equation  (2.2),  the  Addition  Theorem  for  Legendre  Polynomials,  and  Lemma  2.3, 
we  now  obtain 


=  E  E 

n=0  m=-n 


Fn  m(Q!>/^)  ym/A  x\  n 


=  E(E  A«  ^ 

n=0  m——n 

E<  ^Qd"'m(l)  Ynm(0,0)  r"). 

m=l 

Case  I  :  m'  =  0.  Due  to  Lemma  2.3  and  (2.18),  we  have 

- A«<  (^7 J 

=  E(  E  (-i)n'4,  A-  aimJa^'-w  (A)  .y*(m)  r"+ 

E(-i)"'A°n,  •  A-  •  a?Qdg-"'-"*  Q)  Vnm(^,^>)  r"  ) 


n= 0 m=-n 


P  An+n' 


where  the  last  equality  is  obtained  from  Lemma  2.3. 

Case  II  :  m!  <  0.  Using  Lemmas  2.3  and  2.4,  we  have 


r/n'  +  l  -pUzp  yrl  J 

=  E(  E  (-!)"' A?'  ■  A”  aK'aW^n'-l"»l-lm'l  /'I')  .  Ynm(d,<t>)  rn  + 

n*<0  m«-n  V/7/ 

^(-l)n'A^'  •  A™  •  aMam^an+n'-|m'|-m  F„m(^,o)  r"  ) 

m=  1 

X  n  /  *  m'  .  «m  y-m'-m/  ;?v\ 

V'  V'  I  V  ‘*n  »n+n'  la^M  ,  „ 

"2-2-  - -n+n'+l  .  t m  — m' -  K  r  , 

n=0  m=  — n  \  "  ^n+n'  J 


where 


,m'  _  j  (-1)"  ■  if  m  >  0: 

m  \  (-!)"'(- l)m,n<lm'l-lml),  if  m  <  0. 
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Case  III  :  m!  >  0.  From  Lemma  2.3, 


K- 


t 


P 


yi  n'  +  i  -An'  d+pdzp  y^J 

=  E(V  (-1  r'AV'  -A?  (i)  •  rn+ 

n=0  m=-n  Vr/ 

E(-1)”'^'  ■  K  d™'Qd™Qdnz+n'-m-m'  (i)  •  r" ) 

m=  1 


m- 

oo  n 


=  EL  lJ"  t5gMI  CM  A 

n=Om=-n  \ 


„n+n'  +  l  .  j m-m' 
"  '^n+n' 


where 


‘  l  (-!)"  (- 


if  m  <  0; 
jf  m  >  Q. 


As  before,  for  the  terms  with  m  >  0,  we  have  used  Lemma  2.4  to  annihilate  whichever  of  the 
operators  d-  and  d+  occurs  less  frequently. 


The  last  addition  theorem  describes  a  formula  for  expanding  a  spherical  harmonic  of  nonneg¬ 
ative  degree  about  a  shifted  origin.  Its  proof  is  similar  to  those  of  the  first  two  addition  theorems. 
A  more  involved  proof,  based  on  group  representation  theory,  can  be  found  in  [15]. 


Theorem  2.5.  (Third  Addition  Theorem J  Let  Q  =  (p,a,/3)  be  the  center  of  expansion  of  an 
arbitrary  spherical  harmonic  of  nonnegative  degree.  Let  the  point  P  =  (r,9,tp)  with  P-Q  — 
(r1 ,9' ,</>').  Then 


rff'i*,*)  f  ^  JZL-a?  K:;nm  pn  Ynm(°,0) 

___  2^  2-/  TiT - *  n'  —  n  \V-V)r 


n=0  m=  —  n 


AS 


where 


\  -  l)n(  —  l)m,  if  m  m'  <  0; 

J™m  =  ^  (— 1)"(—  if  m  m'  >  0  and  |m'|  <  |m|; 

(-1)",  otherwise. 


(2.19) 


We  are  in  a  position  now  to  develop  the  translation  operators  which  will  allow  us  to  manipulate 
multipole  expansions  and  local  expansions  in  the  manner  required  by  the  fast  algorithm. 

Theorem  2.6.  (Translation  of  a  Multipole  Expansion,)  Suppose  that  m  charges  of  strengths 
.?m  are  located  inside  the  sphere  D  of  radius  a  with  center  a t  Q  =  {p.a.J).  and  that 
for  points  P  =  ( r.9.<p )  outside  D,  the  potential  due  to  these  charges  is  given  by  the  multipole 
expansion 


®(r>  =  £  £  -St  r„”(«U'). 


(2.20) 


n=0  m=-n 


where  P  -  Q  =  (r'.9\  Then  for  any  point  P  =  (r,  9.  d>)  outside  the  sphere  D]  of  radius  (a  ±  p) 


;=  o  k=-} 


(2  21) 


10 


(vlvL'V 


_1i 

.  N 


» 


.  j 


where 


-  E  E 


Ok~m 

J-n 


■  pn  ■  yn~m(^/?) 


AJ 


(2.22) 


with  J*  and  A *  defined  by  equations  (2.14)  and  (2.6),  respectively.  Furthermore,  for  any  p  >  1, 


P  '  Mk  . 

*(P)~E  £  3TT  ■*?(*,*) 

y=o  k——j 


<-  (m)  m* 


(2.23) 


Proof.  The  coefficients  of  the  shifted  expansion  (2.21)  are  obtained  by  applying  the  First  Addition 
Theorem  to  each  of  the  terms  in  the  original  expansion  (2.20).  For  the  error  bound  (2.23).  observe 
that  the  terms  M™  are  the  coefficients  of  the  (unique)  multipole  expansion  about  the  origin  of 
those  charges  contained  in  the  sphere  D,  and  Theorem  2.2  applies  with  a  replaced  by  a  +  p. 

I 

The  second  translation  procedure  is  used  to  convert  a  multipole  expansion  of  the  field  induced 
by  a  collection  of  charges  into  a  local  expansion  inside  some  region  of  analyticity. 

Theorem  2.7.  (Conversion  of  a  Multipole  Expansion  into  a  Local  Expansion,)  Suppose  that  m 
charges  of  strengths  qi,q2,-  ■  ■  ,qm  are  located  inside  the  sphere  Dq  of  radius  a  with  center  at 
Q  =  (p,a,f3),  and  that  p  >  (c  +  l)a  with  c  >  1  (Figure  3).  Then  the  corresponding  multipole 
expansion  (2.20)  converges  inside  the  sphere  Do  of  radius  a  centered  at  the  origin.  Inside  Do,  the 
potential  due  to  the  charges  qi,qv,  -,qm  is  described  by  a  local  expansion: 


J-0 k=-j 


where 


=  g  £  °»  J?  A»  ArYr+n  (a'^ 


Am-k  .  pi+n+1 


(2.24) 


(2.25) 


n= 0  m=-r» 

with  Jf  and  A*  defined  by  equations  (2.17)  and  (2.6),  respectively.  Furthermore,  for  any  p  >  1, 


Lj  ‘  Yk(9,  <(>)  ■  r,+1 


j=o k=-j 


<-  (“)  or 


(2.26) 


Proof.  We  obtain  the  coefficients  of  the  local  expansion  (2.24)  by  applying  the  Second  Addition 
Theorem  to  each  of  the  terms  in  the  multipole  expansion  (2.20).  The  bound  (2.26)  is  an  immediate 
consequence  of  the  simpler  error  bound  (2.4)  and  the  triangle  inequality. 

I 

The  following  theorem  yields  a  procedure  for  shifting  the  origin  of  a  truncated  local  expansion. 
It  is  an  exact  translation,  and  no  error  bound  is  needed. 

Theorem  2.8.  (Translation  of  a  Local  Expansion) 

Let  Q  =  (p,  a,0)  be  the  origin  of  a  local  expansion 

*(P)  =  E  E  O™  '  Y™{0' ,4>)  ■  r'n  ,  (2.27) 

n=0  m=-n 
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3.  The  Fast  Multipole  Algorithm 

In  this  section,  we  present  an  algorithm  for  the  rapid  evaluation  of  the  potentials  and/or 
electrostatic  fields  due  to  distributions  of  charges.  The  central  strategy  used  is  that  of  clustering 
particles  at  various  spatial  lengths  and  computing  interactions  with  other  clusters  which  are  suf¬ 
ficiently  far  away  by  means  of  multipole  expansions.  Interactions  with  particles  which  are  nearby 
are  handled  directly. 

To  be  more  specific,  let  us  consider  the  geometry  of  the  computational  box,  depicted  in  Figure 

4.  It  is  a  cube  with  sides  of  length  one,  centered  about  the  origin  of  the  coordinate  system,  and  is 
assumed  to  contain  all  N  particles  of  the  system  under  consideration. 

Fixing  a  precision  e,  we  choose  p  as  /o^(e)  and  specify  that  no  interactions  be  computed  for 
clusters  of  particles  which  are  not  well-separated.  This  is  precisely  the  condition  needed  for  the 
error  bounds  (2. II), (2. 23)  and  (2.26)  to  apply  with  c  =  2,  the  truncation  error  to  be  bounded  by 
2~p,  and  the  desired  precision  to  be  achieved.  In  order  to  impose  such  a  condition,  we  introduce 
a  hierarchy  of  meshes  which  refine  the  computational  box  into  smaller  and  smaller  regions  (Figure 
4).  Mesh  level  0  is  equivalent  to  the  entire  box,  while  mesh  level  /  +  1  is  obtained  from  level  l  by 
subdivision  of  each  region  into  eight  equal  parts.  The  number  of  distinct  boxes  at  mesh  level  l  is 
equal  to  8*.  A  tree  structure  is  imposed  on  this  mesh  hierarchy,  so  that  if  ibox  is  a  fixed  box  at 
level  /,  the  eight  boxes  at  level  l  +  1  obtained  by  subdivision  of  ibox  are  considered  its  children. 

Other  notation  used  in  the  description  of  the  algorithm  includes 

nearest  neighbor  For  box  i  at  level  /,  a  nearest  neighbor  is  a  box  at  the  same  level  of 

refinement  which  shares  a  boundary  point  with  box  i. 

second  nearest  neighbor  For  box  i  at  level  l ,  a  second  nearest  neighbor  is  a  box  at  the  same 

level  of  refinement  which  shares  a  boundary  po;nt  with  a  nearest 
neighbor  of  box  i. 

$/,,  the  pth- order  multipole  expansion  (about  the  box  center)  of  the  po¬ 

tential  field  created  by  the  particles  contained  inside  box  i  at  level 

/, 

the  p^-order  local  expansion  about  the  center  of  box  i  at  level  /, 
describing  the  potential  field  due  to  all  particles  outside  the  box  and 
its  nearest  and  second  nearest  neighbors. 


the  p,fc-order  local  expansion  about  the  center  of  box  i  at  level  Z, 
describing  the  potential  field  due  to  all  particles  outside  i’s  parent 
box  and  the  parent  box's  nearest  and  second  nearest  neighbors. 


Interaction  list:  for  box  i  at  level  Z,  it  is  the  set  of  boxes  which  are  children  of  the 

nearest  and  second  nearest  neighbors  of  i’s  parent  and  which  are  not 
nearest  or  second  nearest  neighbors  of  box  i. 

Suppose  now  that  at  level  l  -  1,  the  local  expansion  has  been  obtained  for  all  boxes. 

Then,  by  using  Theorem  2.6  to  shift  (for  all  i)  the  expansion  to  each  of  box  i’s  children,  we 

have,  for  each  box  j  at  level  Z,  a  local  representation  of  the  potential  due  to  all  particles  outside  of  j's 
parent’s  neighbors,  namely  Vfij.  The  interaction  list  is,  therefore,  precisely  that  set  of  boxes  whose 
contribution  to  the  potential  must  be  added  to  Vtj  in  order  to  create  Vij.  This  is  done  by  using 
Theorem  2.7  to  convert  the  multipole  expansions  of  these  interaction  boxes  to  local  expansions 
about  the  current  box  center  and  adding  them  to  the  expansion  obtained  from  the  parent.  Note 
also  that  with  free-space  boundary  conditions,  ^o.i  and  'ti,,  are  equal  to  zero  since  there  are  no 
well-separated  boxes  to  consider,  and  we  can  begin  forming  local  expansions  at  level  2. 

Following  is  a  formal  description  of  the  algorithm. 


Initialization 


Algorithm 


Choose  a  level  of  refinement  n  log8  iV,  a  precision  e,  and  set  p  %  log2(e). 

Upward  Pass 

Step  1 

Comment  [  Form  multipole  expansions  of  potential  field  due  to  particles 
in  each  box  about  the  box  center  at  the  finest  mesh  level.] 


do  ibox  =  1, ...,  8" 

Form  a  p'A-degree  multipole  expansion  $n,iboz,  by  using  Theorem  2.2. 

enddo 


Step  2 

Comment  [  Form  multipole  expansions  about  the  centers  of  all  boxes 

at  all  coarser  mesh  levels,  each  expansion  representing  the  potential 
field  due  to  all  particles  contained  in  one  box.  ] 

do  Z  =  n  -  1,  ...,0 
do  ibox  =  1,  ...,8* 

Form  a  ptA-degree  multipole  expansion  by  using 

Theorem  2.6  to  shift  the  center  of  each  child  box’s  expansion 
to  the  current  box  center  and  adding  them  together. 

enddo 

enddo 


Downward  Pass 
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Comment  [  In  the  downward  pass,  interactions  are  consistently  computed 
at  the  coarsest  possible  level.  For  a  given  box,  this  is  accomplished 
by  including  interactions  with  those  boxes  which  are  well-separated 
and  whose  interactions  have  not  been  accounted  for  at  the  parent’s 
level.  ] 

Step  3 

Comment  [  Form  a  local  expansion  about  the  center  of  each  box  at  each  mesh  level  l  <  n- 
This  local  expansion  describes  the  field  due  to  all  particles  in  the 
system  that  are  not  contained  in  the  current  box  or  its  nearest  neighbors.  Once 
the  local  expansion  is  obtained  for  a  given  box,  it  is  shifted,  in  the  second 
inner  loop  to  the  centers  of  the  box’s  children,  forming  the  initial  expansion 
for  the  boxes  at  the  next  level.  ] 

Set  =  *li2  =  •  •  •  =  *1,8  =  (0, 0, ...,  0) 

do  /  =  1, ...,  n  —  1 
do  ibox  =  1,  ...,8* 

Form  *j, n,oz  by  using  Theorem  2.7  to  convert  the  multipole  expansion 
<$i  j  of  each  box  j  in  interaction  list  of  box  ibox 
to  a  local  expansion  about  the  center  of  box  ibox,  adding  these 
local  expansions  together,  and  adding  the  result  to  */,,4OI. 
enddo 

do  ibox  =  1,  ...,8/ 

Form  the  expansion  *j+i,y  for  i&ox’s  children 

by  using  Theorem  2.8  to  expand  'itijhoz  about  the  children’s  box  centers. 

enddo 

enddo 


Step  4 

Comment  [  Compute  interactions  at  finest  mesh  level  ] 
do  ibox  =  1,...,8" 

Form  *j,,6or  by  using  Theorem  2.7  to  convert  the  multipole  expansion 
of  each  box  j  in  interaction  list  of  box  ibox 
to  a  local  expansion  about  the  center  of  box  ibox,  adding  these 
local  expansions  together,  and  adding  the  result  to  * . 

enddo 

Comment  [  Local  expansions  at  finest  mesh  level  are  now  available. 

They  can  be  used  to  generate  the  potential  or  force  due  to  all 
particles  outside  the  nearest  neighbor  boxes  at  finest  mesh  level.  ] 

Step  5 

Comment  (  Evaluate  local  expansions  at  particle  positions.  ] 
do  ibox  =  1, ...,  8n 

For  every  particle  py  located  at  the  point  P}  in  box  ibox, 
evaluate  $n,n,0z{P,)- 

enddo 
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Comment  [  Compute  potential  (or  force)  due  to  nearest  neighbors  directly.  ] 
do  ibox  =  l,...,8n 

For  every  particle  pj  in  box  ibox ,  compute  interactions  with 
all  other  particles  within  the  box  and  its  nearest  and  second 
nearest  neighbors. 

enddo 

Step  7 

do  ibox  =  1, ...,  8n 

For  every  particle  in  box  ibox,  add  direct  and  far-field  terms  together. 

enddo 


Remark:  Each  local  expansion  is  described  by  its  p2  coefficients.  Direct  evaluation  of  this  expan¬ 
sion  at  a  point  yields  the  potential.  But  the  force  can  be  obtained  from  the  gradient  of  the  local 
expansion,  and  these  partial  derivatives  are  available  analytically.  There  is  no  need  for  numerical 
differentiation.  Furthermore,  since  the  components  of  V$  are  themselves  harmonic,  there  exist 
error  bounds  for  the  force  of  exactly  the  same  form  as  (2. 11), (2.23)  and  (2.26). 


A  brief  analysis  of  the  algorithmic  complexity  is  given  below. 


Step  Number 

Operation  Count 

Explanation 

Step  1 

order  Np 2 

each  particle  contributes  to 
one  expansion  at  the  finest 
level. 

Step  2 

order  Np4 

At  the  Ith  level,  8l 
shifts  involving  order  p4 
work  per  shift  must  be 
performed. 

Step  3 

order  <  876 Np4 

There  are  at  most  875  entries 
in  the  interaction  list  for 
each  box  at  each  level. 

An  extra  order  Np 4  work 
is  required  for  the  second  loop. 

Step  4 

order  <  875 Np4 

Again,  there  are  at  most  875 
entries  in  the  interaction 
list  for  each  box  and 
&  N  boxes. 

Step  5 

order  Np2 

One  p<A-degree  expansion  is 
evaluated  for  each  particle. 

Step  6 

order  yiVfcn 

Let  k„  be  a  bound  on  the 
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number  of  particles  per  box 
at  the  finest  mesh  level. 
Interactions  must  be 
computed  within  the  box 
and  its  eight  nearest  neighbors, 
but  using  Newton’s  third  law, 
we  need  only  compute  half 
of  the  pairwise  interactions. 


Step  7 


order  N 


Adding  two  terms  for 
each  particle. 


The  estimate  for  the  running  time  is  therefore 


N  •  (a  ■  (log2(e))2  +  b  ■  (log2(e))4  +  d  ■  k„  +  e)  , 


with  the  constants  a,b,c,d,  and  e  determined  by  the  computer  system,  language,  implementation. 


Remark: 

In  addition  to  the  asymptotic  time  complexity,  asymptotic  storage  requirements  are  an  impor¬ 
tant  characteristic  of  a  numerical  procedure.  The  algorithm  requires  that  $jj  and  'f/ij  be  stored, 
as  well  as  the  locations  of  the  particles,  their  charges,  and  the  results  of  the  calculations  (the  po¬ 
tentials  and/or  electric  fields).  Since  every  box  at  every  level  has  a  pair  of  degree  expansions, 
$  and  ¥,  associated  with  it,  and  the  lengths  of  all  other  storage  arrays  are  proportional  to  N,  it  is 
easy  to  see  that  the  asymptotic  storage  requirements  of  the  algorithm  are  of  the  form 


(a  +  0  ■  p2)  N  ,  or 


(a  +  0  ■  {log2(t))2  N  , 


with  the  coefficients  a  and  0  determined,  sis  above,  by  the  computer  system,  language,  implemen¬ 
tation,  etc. 


Remark:  It  is  clear  that  the  operation  count  for  Step  6  that  the  algorithm  assumes  a  reasonably 
homogeneous  distribution  of  particles.  If  the  distribution  were  highly  non-homogeneous,  then  we 
would  need  to  refine  only  those  portions  of  space  where  the  number  of  particles  is  large.  Although  its 
description  is  more  involved,  an  adaptive  version  retains  both  the  accuracy  and  the  computational 
speed  of  the  algorithm  (see  [4j). 


4.  Conclusions 


An  order  O(N)  algorithm  has  been  constructed  for  the  evaluation  of  Coulombic  fields  and 
potentials  in  large-scale  ensembles  of  particles  in  three  dimensions.  The  algorithm  generalizes  our 
earlier  two-dimensional  results  (see  (7j ,(S])  and  retains  both  the  general  structure  and  the  asymptotic 
complexity  of  its  predecessor.  Here,  we  present  a  description  of  the  procedure  and  its  complexity 
analysis.  Numerical  experiments  are  currently  being  conducted  with  a  computer  implementation 
of  the  algorithm,  and  the  results  of  these  experiments  will  be  reported  in  the  near  future. 


mm*® 


•  •.**  %'  •  ‘  **•*••**  On*  •  Is*  «  .  *  . 


References 

[1]  C.  R.  Anderson,  A  Method  of  Local  Corrections  for  Computing  the  Velocity  Field  Due  to  a 

Distribution  of  Vortex  Blobs,  J.  Comput.  Phys.,  62  (1986),  pp.  111-123. 

[2]  A.  W.  Appel,  An  Efficient  Program  for  Many-body  Simulation,  Siam.  J.  Sci.  Stat.  Comput., 

6(1985),  pp.  85-103. 

[3]  B.  C.  Carlson  and  G.  S.  Rushbrooke,  On  the  Expansion  of  a  Coulomb  Potential  In  Spherical 

Harmonics,  Proc.  Cambridge  Phil.  Soc.,  46(1950),  pp.  626-633. 

[4]  J.  Carrier,  L.  Greengard,  and  V.  Rokhlin,  A  Fast  Adaptive  Multipole  Algorithm  for  Particle 

Simulations,  Technical  Report  496,  Yale  Computer  Science  Department,  1986. 

[5]  A.  J.  Chorin,  Numerical  Study  of  Slightly  Viscous  Flow,  J.  Fluid.  Mech.,  57  (1973),  pp.  785-796. 

[6]  \1.  Danos  and  L.  C.  Maximon,  Multipole  Matrix  Elements  of  the  Translation  Operator,  J. 

Math.  Phys.,  6(1965),  pp.  766-778. 

[7]  L.  Greengard  and  V.  Rokhlin,  A  Fast  Algorithm  for  Particle  Simulations.  Technical  Report 

459,  Yale  Computer  Science  Department,  1986. 

[8]  - ,  A  Fast  Algorithm  for  Particle  Simulations,  J.  Comput.  Phys.,  (to  appear). 

[9]  E.  W.  Hobson,  The  Theory  of  Spherical  and  Ellipsoidal  Harmonics,  Chelsea,  New  York,  1955. 
[10]  R.  W.  Hockney  and  J.  W.  Eastwood,  Computer  Simulation  Using  Particles,  McGraw-Hill, 

New  York,  1981. 

[  1 1]  O.  D.  Kellogg,  Foundations  of  Potential  Theory,  Dover,  New  York,  1953. 

[12]  J.  C.  Maxwell,  Treatise  on  Electricity  and  Magnetism,  3rd  edition  (1891),  Dover,  New  York. 

1954. 

[13]  L.  Greengard,  The  Rapid  Evaluation  of  Potential  Fields  in  Particle  Systems,  Ph.D.  Thesis, 

Dept,  of  Computer  Science,  Yale  Univ.,  1987. 

[14]  V.  Rokhlin,  Rapid  Solution  of  Integral  Equations  of  Classical  Potential  Theory,  J.  Comput. 

Phys.,  60(1985),  pp.  187-207. 

[15]  M.  E.  Rose,  The  Electrostatic  Interaction  of  Two  Arbitrary  Charge  Distributions,  J.  Math,  and 

Phys.,  37  (1958),  pp.  215-222. 

[  16]  P.  R.  Wallace,  Mathematical  Analysis  of  Physical  Problems.  Dover,  New  York,  1984. 


P  —  (r,9,  <t>) 


Figure  1:  Points  P  and  Q  separated  by  a  distance  r\  and 
subtending  an  angle  7  between  them. 


Figure  3:  Source  charges  <71,^2,...,?*  are  contained  in  the 
sphere  D\.  The  corresponding  multipole  expansion  about  Q 
converges  inside  Z?2- 


